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ABSTRACT 


The  work  reported  here  has  been  directed  a long  three  main  avenues. 
First,  a new  method  has  been  developed  for  determination  of  natural 
frequencies  and  associated  mode  shapes  of  thin  elastic  plates  and  also 
shallow  shells  undergoing  free,  undamped  vibration.  This  method  permits 
consideration  of  arbitrary  boundary  conditions  along  each  side  of  a 
polygonal  plate  or  polygonal  plan-form  shell  and  further  offers  a com- 
prehensive evaluation  of  how  well  the  boundary  conditions  have  been 
satisfied  alonq  each  edge.  Second,  the  structural  reliability  of  geo- 
metrically nonlinear  structures  subject  to  stationary  random  excitation 
has  been  investigated  by  two  approaches  which  have  been  found  to  yield 
satisfactory  predictions  of  reliability.  Third,  a technique  has  been 
developed  for  prediction  of  the  statistical  characteristics  of  the 
response  of  a nonlinear  system  to  nonstationary  random  excitation. 
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CHAPTER  I 


VIBRATIONS  OF  THIN  ELASTIC  PLATES 


Introduction 

Problems  involving  vibrations  of  thin,  elastic  plates  occur  in  a 
wide  variety  of  applications  in  structural  mechanics  in  all  aspects  of 
engineering.  An  excellent  comprehensive  summary  of  existing  analytical 
and  experimental  information  pertinent  to  this  problem  area  has  recently 
been  presented  by  A.  W.  Leissa  []].* 

The  present  investigation,  concerned  with  free  vibrations  of  thin 
elastic  plates  is  based  upon  the  assumptions  that  (a)  the  thickness  of 
the  plate  is  small  compared  to  its  other  dimensions,  (b)  there  is  no 
middle  surface  extension  or  contraction,  (c)  points  situated  on  a line 
normal  to  the  undeformed  middle  surface  remain  on  this  same  line  as  the 
plate  deforms  during  vibration,  (d)  stresses  normal  to  the  middle 
surface  of  the  plate  are  neglected,  and  (e)  Hooke’s  law  is  valid. 

Many  problems  involving  free  vibrations  of  plates  have  been  solved 
by  exact  methods,  energy  techniques,  and  various  numerical  approaches, 
in  particular  finite  elements.  The  objective  of  the  present  work  is  to 
present  a more  generalized  approach  that  is  applicable  to  arbitrary  plate 
contours  and  arbitrary  boundary  conditions.  Also,  it  is  desirable  to 
present  a method  that  decreases  the  computer  effort  involved  in  other 
approaches.  With  this  in  mind,  the  present  work  develops  a dynamic 
analog  of  the  Edge-Function  method  originally  due  to  P.  M.  Quinlan  [2] 
and  applied  successfully  by  him  to  a variety  of  elastostatic  problems 
[3], [4]  including  static  deflections  of  thin  elastic  plates  [5]  and  [6]. 

*Numbers  in  brackets  refer  to  references  at  the  end  of  the  text. 


[2] 


Governing  Relationships 

The  fundamental  differential  equation  of  mqtipn  for  fhe  transverse 
displacement  w of  a plate  is: 

D[—  +2  iV_+  l!w]  -t-pjfw  =0  (i„i) 

3x4  3x23y2  3y4  3t2 

where  the  deflection  w (x,y,t)  is  a function  of  the  rectangular  coordinates 
x,  y of  the  plate  middle  surface  and  the  time  t,  anc)  D is  the  flexural 
rigidity  of  the  plate  defined  by: 


12  (1-v2) 

where  E is  Young’s  modulus,  h is  the  plate  thickpess,  v In  Poisson's 
ratio,  and  p is  mass  density  per  unit  middle  surfacq  area  of  the  plate. 
Equation  (1-1)  holds  only  for  small  lateral  motions  of  the  plate  (less 
than  approximately  half  the  plate  thickness)  and  neglects  transverse 
shear  and  rotatory  inertia. 

If  we  introduce  the  Uplacian  Operator  V2  ■jnto  Equation  (1-1),  it 
will  simplify  to: 

DV4w+p— =0  (1-2) 

3t2 

Here  V4  = V2  t V2  and  V2  = — + — 

3x2  3y2 

Using  the  method  of  separation  of  variables  to  solve  Equation  (1-2),  we 
assume  that  vf(x,y, t)  = G(t)  • F(x,y)  (1-3) 

If  we  substitute  Equation  (1-3)  into  Equation  (I*-2),  we  have: 

J G(t)  V4F(x.y)  + F(x,y)  = o 

p 3t2 


(1-4) 


[3] 


Then,  dividing  Equation  (It4)  through  by  G(t)  • F(x,y)  one  finds: 


D v;F(.x,y) , ■ i!Mi)  /G(t) 

p F ( x ,y ) ^4.2 


(1-5) 


The  left  member  of  Equation  (1-5)  is  clearly  independent  of  t,  and 
the  right-hand  side  is  independent  of  x and  y.  Thus,  each  side  of 
Equation  (1-5)  must  be  a constant,  say  K7,  and  we  can  write: 


D vMx^)  _ d2G(t)/dt2  2 

P F(x,y)  ~ G(t) 


(1-6) 


Thus,  the  original  partial  differential  equation  has  been  reduced  to 
two  equations 

V4F(x,y)  - §K2  F{x,y)  = 0 (1-7) 

and 

diGjl).  + K2fi(t)  _ 0 (1-8) 

dt2 

The  solution  of  the  ordinary  differential  Equation  (1-8)  is:. 

G(t)  = A cos(Kt)  + B sin  (Kt)  (1-9) 

where  A and  B~  are  constants,  and  can  be  determined  from  initial  conditions 
To  solve  Equation  (1-7),  one  uses  the  method  of  separation  of 
variables  again.  Assume: 

F(x,.y)  = U(y)  • sin(mx  + a)  (1-10) 

Substitution  of  Equation  ( 1-10)  into  Equation  (1-7)  yields 

dHJk!  . 2m2^UM  + (m-  . K2fr)U(y)  = 0 

dy4  dy2  D 


(1-11) 


[4] 


which  is  an  ordinary  diffqrfntial  equation.  The  sqlution  of  (1-11)  leads 
to  three  different  cases,  a?  follows: 


(1-12) 


(1-13) 

L = [/K2£  - m2  ] 2 

Thus,  the  time- dependent,  fjree,  undamped  vibrations  of  a thin,  elastic 
plate  are  characterized  by  the  lateral  response  w(x,y,t)  given  by  (1-3) 
where  G(t)  is  given  by  (1*9)  and  F(x,y)  is  found  from  ( I - 1 0)  with  U(y) 
given  as  one  or  the  other  of  the  three  cases  discussed. 

The  Edge-Function  method  involves  the  association  of  a different 
system  of  coordinate  axes  with  each  edge  of  the  plate*  Let  us  consider 
a rectangular  coordinate  system  x,y  as  shown  in  Figure  1.  With  respect 
to  this  system,  let  the  vertices  Pl9  P2,  ...  P..  of  the  polygonal  plate 
have  coordinates  (xi,  yi),  (x2,  y2)™ (x^,  y. )— respectively.  Each 
point  P.  is  chosen  as  the  origin  of  a rectangular  coordinate  system 
(x.,  y^),  the  x;.  axis  of  which  makes  an  angle  0.  with  the  positive 
direction  of  the  x-axis  as  shown  below.  This  axis  will  ultimately  be 
chosen  to  coincide  with  an  pdge  of  the  plate. 


if  m2  > 


T5 


U(y)  - 


if  m2  = /plf 

D 


cos  (Ly)  if  m2  < /72P 

K D 


where 


E = [m2  - 


/|/2P  ] ? 


['■>1 


Figure  1 


X 


If  a point  P(x,y)  in  the  (x,y)  coordinate  system  has  coordinates 
x\,  y1#  with  respect  to  the  (x\,  yj  ) coordinate  system,  the  relation 
between  (x,y)  and  (x1#  yx)  is: 


x 

y 


= (x  - x,)  cos  0 + (y  - y ) sin  6 
1 i i i 

= -(x  - x ) sin  6 + (y  - y ) cos  6 

ii  ii 


(1-14) 


as  is  evident  from  Figure  2 


Figure  2 


Analogous  geometric  relations  may  be  found  for  the  coordinates  of  point 

P is  some  (x.,  y.)  coordinate  system  in  terms  of  its  coordinates  in  the 
J J 

( x^ , y. ) system.  Likewise,  the  partial  differential  relations  between 
the  (x,y)  and  the  x. , yj  coordinate  systems  can  be  derived  from  the 

J J 

chain  rule.  For  example,  we  have: 


d_ 

3x 


a 

ay 


a a 

cos  — sin  0. 


J 8xj 


j ayj 


(1-15) 


sin  ej^T  + cos  °j-^T 


The  relationships  between  the  various  second  derivatives  are  found  to  be: 


— = cos2  0 . - 2 sin0.  cose.  — - + sin2  0.  3— 

ax2  J ax.  j j ax. ay.  J ay^ 

J J J J 


>>2  >y2  rx2  3 2 

— = sin2  0 . - - 2 sine . cos0.  + cos2  0.  — 

a.y2  J ax^  j j ax^ay.  J ay^ 


= (y?  ■ aTT^  sin  ei  cos  ei  + (cos20i  - sin2e,)  — 

3x3y  3Xj  3yj  J J J J 3i.3y.  (J.16) 

Comparable  relations  exist  for  the  various  necessary  third  and  fourth 
derivatives. 

The  Edge-Function  Method 

Consider  a thin  elastic  plate  having  n straight  sides  of  lengths 
a1 » a2,  — an.  A right-hand  Cartesian  coordinate  system  in  the  plena 
of  the  plate  is  shown  in  Figure  (3).  The  sides  of  the  plate  make  angles 
Oy  0£9  03»“~“9n  with  the  x-axis. 


[7] 


The  equation  of  undamped  motion  of  this  plate  was  given  as 
Equation  (T-l)  with  respect  to  the  (x,y)  coordinate  system: 

a4w  + 2 s4w  + a^w  + £ a^w  = Q j 

ax4  ax2 ay2  ay4  d at2 


Now,  if  the  vertex  P-  of  the  plate  is  taken  as  an  origin  and  side 

J 

P-  P.  + 1 is  taken  to  coincide  with  the  positive  x.-axis  and  a y.  axis 

vj  \J  J J 

is  selected  which  is  normal  to  x.  and  directed  toward  the  interior  of 

J 

the  plate,  then  a new  coordinate  system  (x.,  y.)  is  obtained.  Since 

J J 

the  bi harmonic  operator  is  invariant  under  rotation  and  translation  of 
axes,  one  immediately  has: 

34  + 2 34  + _3^  a + 2 94  + (1-18) 

3x4  3x23y2  3y4  9x!j  9>u3y*.  9y4. 

J J J J 

Using  this  concept.  Equation  (1-17)  can  be  transformed  as  follows: 

?^T  + ^ + -—  = 0 (1-19) 

8x,  aS^ay^  ay.  D at2 

J J J vJ 


After  applying  the  method  of  separation  of  variables  by  letting 
w(x,  y,  t)  = F(x,y)G(t),  a time-independent  differential  equation  is 
obtained  as  follows: 


[8] 


+ 2 ^ K2F(x,y)  = 0 

95J  3^j  3*j  " 


(1-20) 


or 

v“F(x,y)  - § K2F(x,y)  = 0 


If  the  complementary  solutions  of  Equation  (1-19)  are  F(x-,  yj* 

J J 

j = 1,  2,  — , n,  where  F(x.,  y.)  are  related  to  the  (x.,  y.)  coordinate 

J J J J 

system,  then  from  the  property  of  linear  homogeneous  differential 
equations,  it  follows  that  the  superposition  principle  can  be  applied, 
hence  the  sum  of  those  complementary  solutions  is  also  a solution  of 
Equation  (1-17) , i .e. 

F(x,y)  = l F(Xj,  y.)  (1-21) 

j=l  J J 

Since  F(x-,  y.)  is  a solution  of  Equation  (1-17),  and  Equation 

J J 

(1-20)  is  identical  to  Equation  (1-17),  due  to  the  invariance  property, 
it  can  be  said  that  F(x-,  y,)  is  a solution  of  Equation  (1-20)  or 

J J 


y.i) 


+ 2 


S^Xj,  yj)  dkF(x.,  y.) 

3*j3*5  + SyJ 


- •§  yj)  = 0 

(1-22) 


The  solution  of  Equation  (1-22)  is: 

F(V  *j}  \l}  [AjMe"Hj'yj  + Vj(V]  sin  + aj) 


(1-23) 


where 


[9] 


= 


e“E/j 


if  in.  > /„?p 

3 K D 


if  m.  = /Tp 
J K D 


cos  (L/j) 


if  m.  < /72p' 
J K D 


(1-24) 


and 


H,  = 


[n  • - ] 2 ; 

J K p 


ej  ■ fj  - -FT 


(1-25) 


J. 

Lj  = “ mp2  , m.  is  later  found  to  be  given  by  (5.6). 

u J 

The  required  derivatives  of  (1-24)  with  respect  to  the  (x. , y.) 
coordinate  system,  where  the  origin  is  taken  at  the  vertex  P.  of  the 
plate  and  the  x-  axis  coincides  with  the  edge  P..P.+1  are  readily  found 
through  use  of  the  chain  rule. 


Deflection,  Normal  Slope,  Bending  Moment,  Twisting  Moment,  Kirchhoff  Force 

In  discussing  the  free  lateral  vibrations  of  thin  plates,  the 
deflection,  normal  slope,  bending  moment,  twisting  moment  and  Kirchhoff 
force  can  be  found  as  follows: 

(a)  Deflection 

The  deflection  is  defined  as  w(x,y,t),  where  w(x,y,t)  = 

G(t)  * F(x,y)  or 

n 

w(x,y,t)  = G(t)  • 2 F(x . ,y .) 

j-1  J 3 


(1-26) 


[10] 


Substituting  Equation  (1-23)  into  Equation  (1-26),  it  is  seen  that: 

n °°  - 

w(x,y,t)  = G(t)  • Z Z [AjMe'Hjyj  + B,.Mu^(y)]sin(m.x,.  + as) 


j=l  f1=1 


jMj 


where  u.(y.)  is  defined  in  Equation  (1-24). 

J J 


or 


Thus: 


j j r 
(1-27) 


(b)  Normal  slope  with  respect  to  the  x^-axis: 

The  normal  slope  is  defined  as  6_j(x,y,t)  where: 

e (x,y,t)  = 

3*1 


9.(x,.y,t)  = G(t)  • (1-28) 

3yi 


n 00 

6- (x,y ,t)  = G(t)  • ^ [A jM(mjS i n 9i jCOstm^Xj  + ctj)  - 


HjCos  Oij(sin(Yj  + «j)>  ’ e'HJyJ  + Vj(Yj)] 


(c)  Bending  Moment  with  Respect  to  the  x.-axis: 

, J t 

The  bending  moment  is  defined  as  M (x,y,t)  where 

yi 

M (x.y.t)  = D(32^-t)  + v 32w(Y»t)) 
yi  dy.  3x. 


or 


My  (X,y,t)  = -DG(t)  (3-2F(x’y)  + 

yi  8^ 


32F(x,y)  \ 

V -2  * 


Ak\,  e ; ■ 


[11] 


My _(x,y,t)  = + DC ( t)  J.  [AjMlnK(sinMKj  + vcos  0^j)sin( m j x i + a.) 
1 3 1 

+ 2(l-v)H.m.sin9. .cose. .cos(m.x.+a .)-H.(cos20 . .+vsin20. . ) • 

JJ  ^ J ^ J JJJ  J ^ J ^ J 

sin(mjyaj)}  e'Hjy;i  + BjH5j(Xj’pj)]  {I_29) 

where  U-  are  functions  of  m-,  0..,  a.,  E - , x.  and  y..  In  an  analogous 

J J ' J J J J J 

manner  expressions  are  found  for  the  twisting  moment  and  Kirchhoff 
force,  dentoed  by  R.  A more  general  form  is  obtained  if  Q.  is  taken 
to  represent  these  quantities  where 

^il  = wi 5 Qi2  = 6i’  Qi3  = Hyi  ’ ^4  = ^y.  and  Qi5  = R‘ 

Then  one  has 

‘WV  ■ »<*>  j,  j,  f1-30' 


where  and  are  functions  of  m y 0^.,  ay  Xj,  and  y^,  and  are  given 
in  complete  detail  in  [7]. 


Boundary  Conditions: 

(a)  Simply  supported  edge:  If  the  edge  of  a plate  is  simply 

supported,  the  lateral  deflection  along  this  edge  is  zero,  and  there 
is  no  bending  moment  normal  to  this  edge.  Assuming  that  the  simply 
supported  edge  is  the  i-th  side  of  the  plate,  the  boundary  conditions 
there  are: 

(1)  w(x,y,t)  = 0 at  y.j  = 0 

(2)  My  (x,y,t)  = 0 at  y.  = 0 


(1-31) 


[12] 


(b)  Clamped  edge:  If  the  edge  of  a plate  is  clamped,  the  deflection 

along  this  edge  is  zero,  and  the  normal  slope  of  this  edge  is  also  zero. 

If  the  i-th  edge  is  the  clamped  one,  then  the  boundary  conditions  th^re  are 

(1)  w(x,y,t)  = 0 at  y.  = 0 

1 (1-32) 

(2)  e.(x,y,t)  = 0 at  y.  = 0 

(c)  Free  Edge: 


0)  M (x,y,t)  = 0 at  yi  = 0 

yi 

(2)  R(x,y,t)  =0  at  y.  = 0 


(1-33) 


where 


R(x.y.t)  = [ 


33w(x,y,t) 


+ 


(Z-v) 


33W(x,y,t) 

3x^y. 


which  is  sometimes  termed  the  induced  shear.  The  boundary  conditions 
usually  require  that  a suitable  pair  of  the  general  Equations  (1-30) 
be  specified  on  each  side  of  the  plate.  Physically,  this  means  that  two 
of  the  functions  representing  the  deflection,  the  normal  slope,  the 
bending  moment  or  the  induced  shear  be  defined  on  that  edge.  Let  P.. 
be  any  point  (x..,  0)  on  the  i-th  edge  and  (P^  be  one  of  the  functions 
specified  on  that  side.  Accordingly,  any  boundary  condition  along  the 
i-th  edge  is: 

j,  l Vjn<"pi>  * ViM<pi»  (,-31> 


where  s assumes  any  one  of  the  values  1,  2,  3,  4,  or  5.  Dividing  by  G(t) 
we  have 


cmpi> 


n 

Z 

j=l 


Z 

M=1 


PAjHV'jM*Pi' 


* WV1 


The  identity  (1-35)  is  valid  alone)  the  i-th  edge,  and  must,  be  satisfied 
for  all  x.  in  the  range  (0,  a.)-  Expanding  g(x^)  in  a Fourier  sine 
series,  we  have: 

9<V  ■ j,  CMsin(riii>  ■ jj,  Cfl1HVWP1>  * 

U-36) 

If  we  now  multiply  both  sides  of  (1-35)  by  sin  (—  x. ) , where  N is  any 

i 

fixed  positive  integer)  and  integrate  from  0 to  a..  we  are  led  to  a system 

of  an  infinite  number  of  equations  for  the  values  N~1 , 2,  3, ».  These 

equations  must,  of  course,  be  truncated  after  a finite  number  of  harmonics, 
say  M = L terms.  Then  the  total  number  of  unknowns  and  in 
Equation  (1-34)  is  given  by: 

Total  unknowns  = 2 x (number  of  sides  of  plate)  x L 

(1-37) 

For  the  determination  of  these  unknown  constants,  one  requires  2 x 
(number  of  sides  of  plate)  x L simultaneous  equations.  In  the  usual  plate 
problem,  there  are  two  specified  boundary  conditions  along  each  edge,  thus 
the  total  number  of  specified  conditions  is  2 x (number  of  sides  of  plate). 
If  the  number  N is  L on  every  edge,  then  the  total  number  of  equations  is: 

Total  equations  = 2 x (number  of  sides  of  plate)  x L 

(1-38) 

Comparing  Equations  (1-37)  and  (1-38),  one  observes  that  the  number  N 
of  harmonic  equations  must  equal  the  number  of  terms  in  the  summation 
in  M,  i.e.  it  is  required  that  L = L.  Solution  of  this  set  of  equations 
yields  the  plate  natural  frequencies  as  well  as  relative  amplitudes  of 
free  vibration  in  each  mode. 


[14] 


It  should  he  noted  that  when  a plate  has  a free  corner,  or  if 
the  plate  is  skewed,  additional  functions  are  required.  These  are 
termed  Fractional-Edqe-Functions  and  are  found  from  (1-36)  by  taking 
M as  a fractional  number  and  a.  to  be  non-zero.  The  fractional  value 

J 

is  essentially  determined  by  trial -and-error  so  as  to  minimize  boundary 
residuals  so  as  to  best  satisfy  boundary  conditions.  These  Fractional- 
Edge-Functions  are  added  to  the  Edge-Function  series  (1-34)  in  the 
positions  PI  = L+l  and  L+2  where  L is  the  truncation  value. 

Numerical  Results 

The  first  six  natural  frequencies  for  a square  plate  clamped  on 

all  four  sides  as  obtained  by  the  present  method  (using  various  numbers 

of  harmonics  from  L = 2 to  L = 6)  are  listed  in  Table  1.  Also  presented 

there  are  the  results  of  S.  Tomotika  [8],  D.  Young  [9],  and  $.  Iguchi 

[10]  for  the  same  problem.  The  frequency  parameter  is  X = K 2 /p/D 

a 

where  K is  defined  in  (1-6). 


TABLE  1 


Source 

Mode  and  X 

1 

2 

3 

4 

5 

6 

Tomotika  [8] 

35,99 

— 

— 

-- 

— 

Young  [9] 

35.99 

73.41 

108.27 

131.64 

132.25 

165.15 

Iguchi  [10] 

— 

73.40 

108.22 

— 

132.18 

164.99 

L = 2 

35.572 

72.80 

107.08 

-- 

130.28 

167.78 

L = 3 

35.968 

73.19 

107.08  - 

131.55 

131.79 

163.60 

L = 4 

35.970 

73.35 

108.12 

131.55 

131.79 

164.25 

L = 5 

35.983 

73.38 

108.12 

131 .58 

132.16 

164.85 

L = 6 

35.984 

73.39 

108.20 

131.59 

132.16 

164.93 

[15] 


It  is  interesting  to  note  that  the  use  of  L = 6 involves  only  about  20 
seconds  of  running  time  on  a CDC  3600  computer  to  yield  the  first  six 
frequencies  indicated  in  Table  1.  The  values  of  deflections.,  slopes, 
moments,  etc.  at  any  point  of  the  plate  may  now  be  readily  found  from 
the  appropriate  governing  equations.  There  relative  values  are  obtained 
by  assigning  to  a specified  point  in  the  plate  the  parameter  value  unity, 
then  comparing  the  values  at  other  points  to  this  unit  value.  To  give 
an  indication  of  the  accuracy  of  the  frequencies  tabulated  in  Table  1, 
the  R.M.S.  values  of  the  residual  normal  slopes  along  each  of  the  four 
edges  of  the  clamped  square  plate  when  vibrating  in  its  first  mode  are 
indicated  in  Table  2. 


TABLE  2 

Boundary  residual  of  normal  slopes  corresponding  to  the  first  mode 
of  a square  plate  clamped  on  all  sides. 


No.  of 
Harmonics 

Boundary  of  the  Plate 

AB 

BC 

CD 

DA 

2 

0.0914 

0.0914 

0.0914 

3 

0.0131 

0.0131 

4 

0.0131 

0.0131 

5 

0.0033 

0.0033 

0.0033 

6 

0.0032 

0.0032 

0.0032 

The  first  two  natural  frequencies  of  a square  cantilever  plate  as 
illustrated  in  Figure  4 as  determined  by  the  present  technique,  as  well 
as  found  by  D.  Young  [9]  using  an  energy  approach  are  tabulated  in 


Table  3. 


The  R.M.S.  values  of  a support  deflection  and  free  edge  residual  bending 
moment  along  DA  are  indicated  in  Table  4. 


TABLE  4 
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Results  for  other  cases  such  as  square  plates  simply  supported  on 
all  sides,  rectangular  plates  simply  supported  on  all  sides,  square  and 
rectangular  plates  with  two  opposite  edges  simply  supported  and  the  other 
two  edges  clamped,  square  plates  with  three  edges  clamped  and  the  fourth 
edge  simply  supported,  rectangular  plates  with  two  adjacent  edges  simply 
supported  and  the  other  edges  clamped,  parallelogram  plates  simply 
supported  on  all  sides  and  a rectangular  cantilever  plate  are  given 
in  [7]. 

Conclusions 

It  has  been  found  that  the  Edge-Function  method  is  well-suited  to 
determination  of  natural  frequencies  and  associated  inode  shapes  of  thin 
elastic  plates  undergoing  free  vibration.  For  those  cases  where  other 
methods  have  been  used  by  previous  investigators  to  determine  these 
quantities,  the  Edge-Function  method  yielded  results  in  excellent  agree- 
ment with  existing  analytical  results  as  well  as  experimental  values 
when  those  were  available. 

The  use  of  Edge-Functions  for  investigation  of  natural  frequencies 
and  mode  shapes  of  thin  plates  yeilds  an  additional  bonus  not  ordinarily 
found  in  other  approximate  methods,  namely  rapid  determination  of  how 
well  the  specified  boundary  conditions  have  been  satisfied*  This  is 
done  by  digital  evaluation  of  significant  structural  parameters  along 
the  boundary,  say  deflection  and  slope  for  a clamped  edge,  which  should, 
of  course,  vanish  identically,  but  the  small  residuals  corresponding  to 
each  of  these  parameters  are  routinely  printed  out  be  the  computer.  For 
a simply  supported  edge,  the  small  residual  deflections  and  bending  moments 
at  a number  of  points  along  the  edge  are  routinely  determined.  Thus,  only 
one  may  readily  ascertain  exactly  how  well  the  specified  boundary  conditions 
have  been  satisfied  along  all  boundaries  of  the  plate. 


rial 


CHAPTER  II 

VIBRATIONS  OF  THIN  ELASTIC  SHELLS 

Introduction 

The  objective  of  the  present  investigation  is  to  present  a new 
generalized  approximate  approach  to  the  problem  of  the  free  vibrations 
of  thin,  elastic  shallow  shells.  The  technique  is  applicable  to  shell 
structures  having  an  arbitrary  polygonal  contour  in  the  base  plane  and 
arbitrary  boundary  conditions  along  this  contour.  The  method  extends 
the  Edge-Function  technique  to  the  problem  of  shallow  shell  vibrations 
and  permits  rapid  determination  of  natural  frequencies  and  associated 
mode  shapes  with  only  modest  amounts  of  computer  effort.  Further,  the 
method  permits  evaluation  of  the  errors  involved  in  this  approximate 
approach  by  indicating  boundary  residuals  along  each  straight  boundary 
of  the  base  of  the  shell.  The  specific  example  of  a shallow  spherical 
shell  is  investigated  in  detail. 

Fundamental  Equations 

The  basic  differential  equations  employed  in  this  work  to  describe 
the  behavior  of  a thin  elastic  shell  in  the  absence  of  applied  loads 
are  those  derived  in  Kraus  [11].  The  system  of  equations  reduces  to 
two  simultaneous  differential  equations  for  the  normal  deflection  w 
and  the  stress  function  <f>  of  the  form: 

DV^w  + Vj*  = -ph  (H-n 


Ej-  V2V20  - V2w  = 0 


(II-2) 


I 19 1 


'•'here  V"  is  the  Laplace  operator,  k and  k are  curvatures  in  the  x and 

* y 

2 'S  2 

y directions  respectively,  V.  = k + k , E represents  Young’s 

K y 3x2  X 3y2 

modulus,  h the  shell  thickness,  v is  Poisson's  ratio,  p the  shell  density, 
t denotes  time,  and  D = Eh3/12(l-v2) . Normal  and  shear  stress  resultants 
as  well  as  bending  and  twisting  moments  may  be  expressed  in  terms  of 
w and  4>  as  indicated  in  [11].  The  tangential  displacements  u,  v along 
the  x and  y axes  respectively  may  be  expressed  in  terms  of  w and  i|>  from 
relations  presented  in  [4].  For  example,  the  first  of  these  relations  is: 


(11*3) 


Let  us  consider  tfie  case  of  a shallow  spherical  shell  for  which 
kx  = ky  - ky=  k = 1/R  where  R is  the  radius  of  curvature  of  the  shell 
middle  surface.  If  both  w and  d>  exhibit  a harmonic  dependence  on  time, 
i .e.  if 

w = W(x,y)  sinwt  (I 1-4) 

<J>  = 4>(x,y ) sinwt  ( 1 1-5 ) 

then  Equation  (1)  and  (2)  reduce  to 


DV4I/I  + kV2(J)  = phw2W  (1 1-6) 

* kv2w  = o (n-7) 

where  is  the  natural  circular  frequency  of  free  vibration.  It  is 
possible  to  uncouple  ( 1 1-6)  and  (1 1-7)  into  the  form 


v6w  - o4v2w  = o 

Ve<(>  - p“V4(j)  = 0 


( 1 1-8) 
(II-9) 
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where 

p"  = ^(pio2  - Ek2)  (11-10) 

Since  there  are  four  conditions  to  be  specified  on  each  edge  of  the 
polygonal  boundary  the  solutions  for  w and  <j>  must  contain  four  arbitrary 
constants. 


For  convenience,  let  us  set 

\p  = V4W  - p4W 

(11-11) 

Substitution  of  this  relation  into  (8)  leads  to 

V2\p  - 0 

(11-12) 

Formulation  of  Solution 

The  projection  of  the  shell  onto  the  x-y  plane  is  considered  to  be 
a two-dimensional  convex  simply-connected  region  R,  the  closed  boundary 
B of  R being  a polygon  of  J sides.  As  in  Figure  5,  a set  of  rectangular 
Cartesian  coordinate  axes  Oxy  is  chosen  and  the  vertices  and  sides  of 


Figure  5 


the  polygon  B are  numbered  1 to  J.  The  typical  vertex  P.  has  coordinates 
( X . , Y.)  in  the  Oxy  reference  frame  and  the  typical  side,  the  j**1,  is 

J J 

P.P.+l,  having  length  a-  and  making  an  angle  <J>.  with  the  positive 

J J J J 

direction  of  the  Ox  axis.  Associated  with  each  side  j is  a set  of 
edge-axes  x.y . , having  origin  at  P.  and  directed  along  the  (inward) 

JJ  J 

normal  to  that  side. 

The  fundamental  series  solution  of  Equation  (11-12)  in  rectangular 

coordinates,  derived  using  the  separation  of  variables  technique  is: 

tjj  = >:  {[Asin(mx  + a ) + B cos(mx  + 3x)]e”my  + [C  sin(mx  + a ) 

m m in  m m 

+ Dmcos(mx  + £x)]e+my}  (11-13) 

where  m,  Am,  B^,  Cm,  and  Dm  are  as  yet  undetermined  constants.  The 
phase  shifts  am  and  3m*  although  unnecessary  here,  are  included  for 
sake  of  generality  later.  Since  the  Laplacian  operator  V2  is  invariant 
under  transformations  which  produce  translation  and  rotation  of  axes, 
it  follows  that 


= 0 


(11-14) 


The  solution  of  this  equation  will  then  take  the  form 

4>  = l {[A^sin(m.x,  + a J + B|jcos(ni.x.  + ej]e_nijyj 
m j j m m j j m 

J 

+ [C^sin(ri)jXj  + o^)  + D^cos(m^xj.  + Bm)]e+mjyj  (11-15) 

JL  U. 

in  terms  of  the  coordinates  (x.y.)  of  a point  referred  to  the  j system 

J J 

of  axes.  Subscripts  and  superscripts  j have  been  added  to  m and  the 

J.U 

constants  Am>  Bm,  Cm,  and  Dm  to  indicate  association  with  the  j side 


[22] 


of  the  boundary  B.  Clearly,  similar  expressions  arise  for  solutions 
associated  with  each  of  the  J sides  of  the  polygon  and  all  such  solu- 
tions may  be  superposed,  the  governing  differential  equation  being 
linear.  Thus,  a general  solution  may  be  put  in  the  form 

j 

* = ^ (m£  } .{[Amsin(mjxj  + °Sn>  + Bmcos(mjxj  + *V*r'njyj 
J 

+ [C^sin(m.x.  + aj  + D^coslnv-Xj  + 3m)]emjyj)  (11-16) 

Analogous  to  the  solution  for  a simple  rectangular  region,  the  constants 
m-  in  (11-16)  are  chosen  so  that 

J 

m.  - 2Mrr/a.  (11-17) 

J J 

for  M a non-negative  integer,  in  order  to  facilitate  generation  of  the 
boundary  condition  equations  later.  Obviously  the  function  ^ contributes 
to  the  determination  fo  the  displacement  components  and  stress  couples. 
Thus,  in  general,  in  formulating  the  boundary  conditions  for  the  typical, 

Ml  i,L 

j , side,  the  coefficients  in  (11-16)  associated  with  the  j n side 
can  be  made  implicitly  dependent  on  the  relevant  boundary  conditions. 

In  the  general  solution  at  any  point  within  the  region  R,  the  displace- 
ment components  and  stress  resultants  will  be  composed  of  contributions 
from  each  side  of  B,  each  contribution  being  dependent  on  the  boundary 
conditions  imposed  on  that  particular  side.  Thus,  invoking  St.  Venant's 
principle,  it  is  physically  desirable  that  such  contributions  decay  with 
increasing  distance  from  the  corresponding  side.  So,  in  11-16),  the 
coefficients  and  are  set  equal uto  zero.  Consequently,  tJj  is  chosen 


to  be: 


[23] 


J 

ip  = z 

j=l 


sCA^sinCnijXj  + <*„,)  + B^costnKX.  + ^)]e'mJyj  (11-18) 


in  which  m.  is  defined  in  (11-17). 

J 

Each  of  the  fundamental  terms  e“mjy.isin(mixi  + a)  and  e~mJyJcos(m  .x . 

j j j j 

+ ^ in  ijj  now  possesses  the  advantageous  properties  of  being  directly 
associated  with  the  j-th  side  of  the  boundary  and  of  decaying  in  contri- 
bution to  the  overall  solution  with  increasing  distance  y.  into  the 

J 

region  from  that  side.  Such  functions,  similar  both  in  form  and  notation 
to  those  employed  by  Tai  and  Nash  [7],  are  termed  edge-functions. 
Comparable  expressions  for  may  be  written  but  are  omitted  for  brevity. 

If  one  substitutes  (11-18)  into  (11-11)  and  also  carries  out  a comparable 
procedure  for  <{>,  one  finally  obtains 


= z z{r-l*-Aj  e'mJyj  + Aj  e‘(mj  + P2)^] 
„tL  4 AlMe  A3Me  3 J 


j=l  M p4  IM 
»?  . n2\l/2 


A4Me"^J'  ' P ^ ^Isin^.Xj  + o^)  + 

[-^  B’j Me’nl Jyj  + B3He-(lnj  + + 


B4Me’("lj  “ P ) ' yj]cos(mjxj  + Bm)} 

♦ = s mjyiAlM)e"mjyj  + 

j=l  H P4  2m. k J J 1M 


and 


j=l  M p^  tn  2m‘rk 

J 


* p2)V2^  - ^ - p2),/2^ 


Ehk 


•sinOyj  +am)  + [-7(B2M-^Vj 

. e-mJy3  + Ehk(Bj  e-(m!  + P2)1/2yj  . 


B1M) 


p2  v 3H 

„?  . „2,1/2w_. 


L 2 1 /2 

B4Me"(rnj  " P ^ yj)]cos(m.jx.j  + em)} 


(H-19) 


(11-20) 


[24] 


The  expression  for  the  transverse  deflection  (11-19)  reduces  in 
the  case  of  zero  curvature,  k = 0,  to  the  corresponding  expression  obtained 
by  Tai  and  Nash  [7]  for  the  transverse  deflection  of  a thin  flat  plate. 

Shell  boundary  conditions  are  expressed  in  terms  of  various  combinations 
of  the  parameters  W,  <&,  the  in-plane  displacements  uq  and  v^  along  and 
perpendicular  to  the  q-th  side  of ' the  boundary,,  the  rotation  3 about  the 
tangent  to  the  q-th  side,  the  rotation  3^  about  the  normal  to  the  q-th 
side,  the  in-plane  stress  resultant  Nq,  the  bending  moment  about  the 
q-th  side,  the  twisting  moment  about  that  same  side,  the  Kirchhoff 
effective  shearing  stress  resultants  Tq  and  Vq,  and  the  in-plane  shearing 
stress  resultant  all  associated  with  the  q-th  edge.  Expressions 
for  each  of  these  twelve  parameters  are  derivable  in  terms  of  the 
coefficients  and  but  are  omitted  for  brevity.  However,  it  is 
convenient  to  symbolize  these  boundary  functions  in  the  form 

At 

where  t is  an  integer  ranging  from  1 to  12  for  the  above  parameters  W, 

$,  u , etc.  respectively.  Further,  let  us  introduce  the  function 
to  denote  the  edge- function  contribution  toAj.  stemming  from  the  j-th 
coordinate  system  and  associated  with  the  particular  value  M.  In  this  case 
each  of  the  above  boundary  functions  may  be  written  in  the  form 
J i 

At 53  s E AtEM  (n-2i) 

1 j=l  M Z M 

It  is  necessary  to  determine  values  for  the  frequency  as  well  as 
the  coefficients  and  so  as  to  satisfy  boundary  conditions.  Thus, 
boundary  conditions  require  that  four  of  the  functions  be  specified 
on  each  edge  of  the  polygonal  boundary.  This  is  accomplished  as  follows. 
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If  P‘  is  any  point  (x  ,0)  on  the  q-th  side  of  the  boundary,  as  in 
Figure  1,  the  function  ^t(Pq)  i$  specified  for  four  values  of  the 
parameter  t t/j,  t^  and  t^,  say.  Consequently,  equation  (11-21) 

implies  that 


vp;>  ■ jf,  ; 

1 = Eq’  tq’  tq’  tq;  ^ B 1 *2 J 


(11-22) 


In  order  to  obtain  identity  equations  for  the  coefficients  and  B^, 
it  is  noted  that  the  series  of  terms 

(M)  t M q 

in  the  summation  (11-22)  does  not  contain  negative  exponential  factors 
since  y = 0,  so  that  each  term  involves  only  x . Since  the  negative 

H 4 

exponentials  tend  to  diminish  boundary  influences  this  series  of  terms 
is  dominant  and  consequently  (11-22)  may  be  written 

MPq>  ' ^ j AtEM(Pq5  = J { 4f,  (hiAiM  + 9iBiM> 


M i~l 


sin(Vq  + %)  + ^ <hi+rA?M  + gi+4B^) 
cos(mqxq  + em)} 


(11-23) 


The  dash  in  the  summation  2 , indicates  that  the  term  j = q is  omitted. 

j=l 

The  factors  h and  g in  (11-23)  are  independent  of  xq  and  may  readily  be 
obtained  from  the  expressions  for  tE^j  by  setting  yq  = 0.  Thus  the  left 
hand  side  of  (11-23)  is  some  function  H(xq)  of  xq  and  the  identity  can  be 
satisfied  provided  that  H(xq)  can  be  expanded  in  the  trigonometric  form 
of  the  series  on  the  right-hand  side.  On  choosing 


mq  = 2Mir/aq,  M = 0,1,2,... 


(11-24) 
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"in  = ‘in  = 0 


(11-25) 

the  right-hand  side  of  (11-23)  may  be  considered  the  Fourier  series 


expansion  of  H(x  ) in  the  range  [0,a  ]. 

M H 

Several  methods  of  approximately  satisfying  the  boundary  identities 
(11-23)  may  be  employed.  Perhaps  the  most  obvious  is  to  multiply  (11-23) 
by 

2N 

sin  x dx 

aq  * ^ 

and  integrate  from  0 to  a , then  repeat  using  the  cosine  function. 

This  leads  to  a set  of  equations  of  the  form 

G(xq)  = 0 (11-26) 

which  is  actually  an  infinite  set  of  equations  in  an  infinite  number 
of  unknowns.  It  is  of  course  necessary  to  truncate  these  equations  at 
some  level,  say  L-,  for  each  side  j of  the  polygon. 

J 

Because  the  truncation  of  the  harmonic  series  expansion  of  G(xq) 
at  N = Lq  implies  that  the  approximation 


G(V  = K + tENC0S  T*  xq  + FNsin  V 


(11-27) 


is  being  used,  another  way  to  approximate  G(xq)  is  to  use  trigonometric 
interpolation  to  fit  a finite  trigonometric  series  to  G(x  ) by  a discrete 
least  squares  method.  As  is  shown  in  [12],  if  a discrete  Fourier  series 
interpolation,  based  on  discrete  least  squares  fitting  at  2k*-l  equidistant 
internal  points  in  the  interval  [0,aq]  is  used  in  approximation  (11-27), 
then 


2k* 


en  = f k!0  wkG(xk)cos  xk 


Fm  = 


1 


2k*-l 

i 

k= 


N ' k*  wkG<xk)sin  xk 


= 0,1 ,2,. . . ,L 


= 1,2 L, 


(11-28) 
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where 

xk  = kaq/2k* 

and  are  weight  factors  defined  by 

wk  - 1/2:  k = 0,  k = 2k* 

=1:  0 < k < 2k* 


(11-29) 


(II- 30 ) 


The  expressions  (11-28)  for  and  may  be  set  to  zero  thus  generating 
the  requisite  set  of  simultaneous  boundary  equations.  This  interpolation 
method  of  setting  up  the  simultaneous  equations  has  the  advantage,  over 
the  first  method  above,  of  requiring  less  computational  time.  Furthermore, 
it  is  found  in  practice  that  there  is  very  little  difference  in  accuracy 
between  the  results  obtained  from  either  method. 

The  set  of  equations  arising  from  expression  (11-27)  in  which  the 
coefficients  of  the  expansion  must  be  set  to  zero,  ensures  that  the 
boundary  conditions  (11-22)  are  satisfied  at  all  points  of  the  q-th  side, 
q = 1,2,...,J,  except  possibly  at  the  vertices.  If  the  boundary  conditions 
are  not  automatically  satisfied  at  the  end  points  of  each  side  of  the 
polygon,  they  must  be  imposed  at  these  points.  Thus  it  is  required  to 
set  to  zero  both  the  boundary  function  and  its  derivative  along  the  relevant 
side,  so  that 

G(xq)  = 0 (II-31) 

and 

f^G(xq)  = 0 (H-32) 

for  both  ends  xq  = 0 arid  xq  = aq  of  the  interval  and  for  each  of  the 
boundary  functions  corresponding  to  t = tf  t"",  t^,  and  tf  The  effect  of 
these  point  equations,  which  are  referred  to  as  vertex  equations,  is  not  alone 


to  ensure  that  the  boundary  conditions  are  satisfied  at  the  end  points  of 
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the  interval,  but  also  to  increase  the  rate  of  convergence  of  the  Fourier 
series  from  order  1/M  to  at  least  order  1/M  . Such  equations  clearly 
give  rise  to  the  necessity  of  introducing  additional  unknowns  in  the 
basic  solutions  (11-19)  and  (11-20).  Since,  in  general,  there  are  four 
vertex  equations  arising  for  each  boundary  condition  on  each  side,  a total 
of  sixteen  additional  unknowns  associated  with  each  side  are  required. 

These  are  provided  by  using  fractional -edge-functions.  The  latter  are 
edge  functions  of  the  types  in  expressions  (11-19)  and  (11-20)  within 
the  summation  signs.  Whereas  in  the  above  harmonic  edge-functions 

mq  = 2M7r/aq,  M = 0,1, = 0;  $m  = 0 (11-33) 

for  the  fractional -edge- functions  new  values  of  these  parameters  are 
defined,  so  that 

V 2V'V  (II-34) 

where  VM  is  non-integral.  The  sixteen  requisite  unknowns  can  be  provided 
by  including  in  the  solution  two  fractional-edge-functions  generated 
simply  by  choosing  two  different  values  for  V^.  They  may  be  incorporated 
inthegeneral  solution  by  adopting  the  convention  that  they  correspond  to 
M - -1  and  M = -2  in  the  summations  in  (11-19)  and  (11-20).  It  may  be  noted 
that,  in  theory,  the  solution  is  independent  of  the  choice  of  the  vertex 
numbers  and  of  am  and  for  the  fractional -edge- functions.  However, 
in  practice  it  is  found  that  a judicious  choice  of  values  will  increase 
the  convergence  of  the  truncated  series  solution. 

It  was  found  in  [7]  that  the  problem  of  free  vibrations  of  flat  plates 
could  be  successfully  treated  by  use  of  a combination  of  edge  functions  and 
fractional  edge  functions.  That  is,  the  prescribed  boundary  conditions 
could  be  satisfied  to  the  desired  degree  of  accuracy.  In  the  case  of  shell 
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vibrations  this  is  not  the  case  and  it  becomes  necessary  to  employ 
not  onlv  the  two  types  of  functions  used  in  plate  analysis,  but  in 
addition  a third  type  of  function  termed  a shell  polar  function. 
This  function  contributes  to  the  time-dependent  deflections  in  the 
interior  regions  of  the  shell  away  from  the  boundaries,  whereas  the 
edge  functions  contribute  largely  in  the  immediate  vicinity  of  the 
shell  edges.  The  shell  polar  functions  are  of  the  form 

<j>  = Re[EZX  + F ZZX+1] 


(11-35) 


W = RE [F 1 , ZX] 


where  z = x + iy  (origin  at  any  convenient  interior  point),  z is  the 
complex  conjugate,  E,  F,  and  F'  are  complex  constants,  and  X is  an 
integer.  Details  of  the  derivation  and  accompanying  properties  of 
these  special  functions  are  presented  in  the  Appendix. 

The  shell  polar  functions  are  appended  to  the  solutions  to  replace 
some  (or  possibly  all)  of  the  fractional  edge  functions.  The  selection 
of  the  ratio  of  fractional  edge  functions  to  shell  polar  functions  is 
made  on  the  basis  of  experience  and  judgment  in  use  of  this  technique 
but  with  the  requirement  that  boundary  conditions  be  satisfied  to  the 
prescribed  degree  of  accuracy.  In  the  present  investigation  of  a combina- 
tion of  half  edge  functions  and  half  polar  functions  was  found  best  for 
satisfying  the  prescribed  boundary  conditions. 

Application  of  the  boundary  conditions  leads  to  a system  of  homo- 
geneous equations  in  the  unknowns  and  of  the  form 


SCab{u))Ab  = °>  b = 1,  2 


(11-36) 
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where  represents  A:^  and  B^.  The  C ^ are  coefficients  obtained  from 
the  interpolation  method  mentioned  previously  and  are  functions  of  the 
unknown  frequency  u>.  The  frequencies  are  found  by  settinq  the  determinant 
of  the  system  equal  to  zero,  viz: 

|Cab(u)|  = 0 (11-37) 

Values  of  co  are  determined  by  a suitable  iteration  method,  such  as  the 
bisection  procedure. 

Numerical  Results 

Let  us  consider  the  free  vibrations  of  a thin,  elastic  shallow 
spherical  steel  shell  of  square  plan-form.  The  radius  of  curvature  is 

30.0  inches,  the  shell  thickness  0.05  inches,  the  length  of  each  side 

12.0  inches  (in  plan-form),  E = 30  x 106  lb/in^,  Poisson's  ratio  = 0.3 
and  the  material  density  is  0.732  x 10”  lb-sec  /in  . Boundary  conditions 
to  be  considered  are  (a)  all  sides  simply  supported,  and  (b)  all  sides 
clamped.  Natural  frequencies  and  associated  mode  shapes  are  desired. 

The  natural  frequencies  are  indicated  by  Equation  (11-37)  and  a computer 
program  for  determination  of  these  frequencies,  associated  mode  shapes, 
and  boundary  residuals  is  available  upon  request  from  the  authors  [13]. 

For  the  case  of  simply  supported  edges,  the  first  four  natural 
frequencies  are  indicated  below.  Results  are  compared  to  those  due  to 
Vlasov  [14]  as  well  as  Malkina  [15].  The  investigation  [15]  applies  to 
the  case  of  a spherical  dome  with  circular  plan-form  but  is  employed 
here  as  an  approximation  by  considering  the  greatest  circular  base  to 
be  inscribed  in  the  square  plan- form  shell  under  consideration. 
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Frequency  1234 

Present  Method 

L=1  6748.86  6759.30  6763.15  6829.48 

L=7  6748.85  6759.30  6863.63  6795.35 

Vlasov  [14]  6749.51  6759.29  6827.35 

Malkina  [15]  6749.42  6762.31 

The  computer  program  developed  indicates  the  small  residual  bending 
moments  along  each  of  the  sides  of  the  shell.  For  example,  for  the 
first  natural  frequency  if  the  peak  bending  moment  at  the  mid-point  of 
the  shell  is  taken  to  be  unity,  the  root-mean-square-boundary  residual 
bending  moment  (for  L.^7)  is  found  to  be  0.37  x 10”  . This  same  function 
for  the  fourth  natural  frequency  is  found  to  be  0.88  x 10”5.  Clearly 
these  constitute  very  adequate  satisfaction  of  boundary  conditions. 

For  the  case  of  a clamped  edge  spherical  shell,  the  first  four 
natural  frequencies  are  indicated  below,  as  well  as  those  found  by  the 
Malkina  method  [15]  again  usinq  the  greatest  inscribed  circular  base. 

Frequency  1234 

Present  Method 

L=1  6746.24  6794.72 

L=7  6746.24  6788.39  6992.20  7135.09 

Malkina  [15]  6746.24  6796.89  6998.91 

Again,  the  computer  program  displays  the  small  residual  deflections 
along  the  shell  boundary.  For  example,  for  the  first  natural  frequency 
the  root-mean-square  boundary  deflection  is  0199  x 10  (for  L=3)  and 
for  the  fourth  natural  frequency  it  is  0.75  x 10"^  (for  L=3)  compared 
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to  peak  mid-point  deflection  taken  to  be  unity.  Mode  shapes  are  found 
by  calculating  the  relative  deflections  of  closely  spaced  points  on  the 
shell  then  connecting  all  zero-deflection  points  at  a given  frequency. 

Conclusions 

The  present  investigation  has. indi cated  that,  for  the  geometry 
considered,  the  Edge  Function  Method  is  well-suited  to  determination  of 
natural  frequencies  and  associated  mode  shapes  of  thin  elastic  shallow 
shells  undergoing  free  vibration.  For  the  cases  discussed  through 
specific  examples,  the  present  analysis  yielded  results  in  excellent 
agreement  with  existing  analytical  results. 

The  computer  program  developed  during  the  course  of  the  present 
investigation  is  applicable  to  any  spherical  shell  of  n-sided  plan-form. 
Further,  the  boundary  conditions  along  each  edge  are  completely  arbi- 
trary and  may  be  different  on  the  various  edges  without  giving  rise  to 
complications  in  frequency  determination. 

The  use  of  the  present  technique  of  analysis  for  determination  of 
natural  frequencies  'and  mode  shapes  of  thin  elastic  shells  yields  an 
additional  bonus  not  originally  found  in  other  approximate  methods,  namely 
rapid  determination  of  how  well  the  specified  boundary  conditions  have 
been  satisfied.  This  is  accomplished  through  digital  evaluation  of  sig- 
nificant structural  parameters  along  the  boundary,  say  delfection  and 
slope  for  a clamped  edge  which  should,  of  course,  vanish  identically  but 
the  small  residuals  corresponding  to  each  of  these  parameters  are  routinely 
printed  out  by  the  computer  and  root-mean-square  values  along  any  of  the 
polygonal  edges  displayed.  Thus,  one  may  readily  ascertain  exactly  how  well 
the  specified  boundary  conditions  have  been  satisfied  along  all  boundaries 
of  the  plate. 
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APPENDIX 

SHELL  POLAR  FUNCTIONS 


Equation  (1 1-9)  gives  for 

V4(V4-  p)c|)  = 0 (A-l ) 

one  solution  of  which  is 
V4(|>  = 0 

from  which,  following  O'Callaghan  [12] 

<t>  = Re[EzA  + FzzX+1];  z = x + iy  (A-2) 

where  E and  F are  complex  constants  and  A is  arbitrary.  Accordingly  on 
substituting  for  V4cf>  in  (7),  we  obtain 

V2W  = 0 (A-3) 

and  Equation  (II-6)  then  gives 

W = — — v2<t>  (A-4) 

ph(ja2 

On  introducing  the  operators 


_ .id)  9 • 

= ie  q — ~ ie 


9z 


from  which 


q*i 


(A-5) 


■ JL  + 1L=  4 _lL_ 

gx2  gy2  3Z3z 

and  on  operating  with  V2  on  (A-2)  and  substituting  in  (A-4),  we  obtain 

W = 4-Cx-+-U  Re[FzX]  = Re(F'zX] 
phrxo2 


(A-6) 
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where 

F*  = qF;  q = -4-^li  (A-7) 

phm2 

To  fit  into  the  general  polar  form  for  the  t ^ derived  function 

f)t  = Re[EAzX*  + F(BzzX*+1  + CzX*+2)]  (A-8) 

where  A,  B,  and  C are  complex  functions  of  <j>  and  A,  and  where  necessary 


will  have  a subscript  t - B^,  C^.  - to  distinguish  the  different 

functions.  Equation  (A-6)  for  W then  gives: 


A 

II 

CQ 

• B 

o 

II 

0; 

C = 2;  A*  = A- 2 

(A-9) 

Values 

of  A,  B, 

and  C 

for  each  of  the  remaining  eleven  parameters  cj>. 

uq’  vq' 

, etc.  are 

tabulated  below. 

Function  t 

A* 

At 

Bt 

Ct 

w 

1 

A- 2 

q 

0 

0 

<l> 

2 

A 

l 

1 

0 

Uq 

3 

A-l 

Aoe1^ 

Boe^ 

vq 

4 

A-l 

iA3 

iB3 

"lC3 

Bq 

5 

A-l 

A4/r 

B4/r 

C4/r-i  C-je^q/z 

Nq 

6 

A- 2 

(X-l)e2iA| 

(X+l)e2i4>q 

2(X+1) 

V 

7 

A-4 

0 

0 

D*A(A-l)e2l<l>q 

8 

A-4 

0 

0 

iC7- 

Tq 

9 

A- 2 

a12 

B12 

C12+cg/rz2 

Vq 

10 

A- 5 

0 

0 

(X-2)i<},qC8/r 

Sq 

11 

A-l 

Vr 

B4/r 

(yr-Ae^q  Z}/2 

Nnt. 

12 

A- 2 

-iX(X-l)e2l<J,q 

-i  XX ' 

0 

D*  = D(l-v) 

X'  = (X+l  )e2icf> 

q 
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CHAPTER  III 

RELIABILITY  OF  NONLINEAR  STRUCTURES  UNDER  RANDOM  EXCITATION 
Introduction 

The  reliability  of  structures  subject  to  random  excitation  has 
become  a topic  of  increasing  importance  to  designers.  For  example,  flight 
vehicles  are  often  subject  to  turbulence,  random  wind  gusts,  jet 
engine  pressures,  rocket  exhausts,  and  other  pressure  fields  that  are 
random  in  nature.  Actually,  these  are  almost  always  nonstationary 
phenomena.  However,  because  of  the  difficulties  involved  in  analysis 
of  nonstationary  situations  it  is  necessary  to  first  consider  the 
response  of  structures  to  stationary  excitation,  i.e.  the  statistical 
characteristics  of  the  pressure  fields  are  invariant  under  time  shifts. 

The  nonlinearity  under  consideration  is  geometric  in  nature,  which 
characterizes  large  deflections  of  structural  members.  The  present 
investigation  is  confined  to  the  response  of  a single  degree  of  freedom 
system. 

In  the  present  work  the  perturbation  method  [16]  is  used  to  deal 
with  the  nonlinear  differential  equations.  Since  this  method  has  its 
limitations  [17],  we  also  use  the  Monte  Carlo  simulation  method  [18] 
to  ascertain  the  range  of  validity  of  the  perturbation  results.  Knowing 
the  response,  it  is  usually  not  difficult  to  determine  significant  system 
stresses.  From  this  type  of  analysis  together  with  considerations  of 
first-crossing  and  fatigue  [19]  it  is  possible  to  predict  the  reliability 
of  narrow-band  type  structures  subject  to  stationary  random  excitation. 

The  Perturbation  Method 

If  a nonlinear  single-degree-freedom  system  is  governed  by  the 
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equation: 

X(t)  + Zuin£X(t)  + u>2(X(t)  + yX3(t) ) = f(t)  (III-l) 

where  c is  the  damping  ratio,  wn  the  natural  frequency  of  the 
corresponding  linear  system,  y a perturbation  parameter  which  is 
dependent  upon  properties  of  structure,  and  f(t)  is  a generalized  force 
which  is  a stationary  Gaussian  random  process  with  mean  zero  and  the 
spectral  density  function 

We  assume  a solution  of  (III-l)  in  the  form  of  a power  series 

in  : 

X(t)  " XQ(t)  + yX-j(t)  + y2X2(t)  + ... 

For  convenience,  let 

X = X0  + X1  + v2X2  + ...  {II 1-2) 

If  E [ ] denotes  ensemble  average,  then 

E(X  X*]  = E[XqX*]  + u(E[X0X|]  + E[X*X-|]) 

+ u2(E[X0X|]  + E[X1X*3  + e[x*x2]  + ...  (Ill -3) 

where  X*  = X(t+t)  and  E[X0X0],  E[XQXf],  EtXJX,],... 
are  to  be  determined. 

When  ( I I I -2 ) is  inserted  in  (III-l),  we  get 

(XQ  + X-j  + y2X2  +...)+  2(unc(X0  + y X*j  + y^X2  + ...) 

+ w2  [X0  + uX1  + p2X2  + ...  + v(X0  + vX1  + ...)3]  = f(t) 
or 

(*o  + 2“n^o  + unV  + + 2(VX1  + + 

+ u2(X2  + 2%d2  + w2X2  + 3X2X-|Co2)  + u3  (...) 

+ ...  = f(t)  (II 1-4) 

The  coefficient  of  each  power  of  v must  separately  vanish  since  (III-4) 
is  to  be  satisfied  identically  in  u.  Thus,  we  obtain  a set  of  governing 
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equations  for  the  X(t).  The  first  of  these  represents  the  unperturbed 
system. 


*o  + 2V^  Xo  + “nXo  = 
h + 2“n«Xl  + “nXl  = '“nXo 

X2  + 2a)ncX2  + u2X2  = -3a12X2X1  (II 1-5 ) 

Because  we  treat  the  excitation  f(t)  as  a stationary  random 
process,  the  response  is  also  a stationary  random  process  after  the 
transient  period,  and  properties  of  the  response  can  be  deduced  from 
that  of  the  excitation.  From  random  process  analysis  [20],  we  obtain 
the  solution  of  (II 1-5 ) and  also 


E[X0X*]  = 
E[X0X*]  = 

E(x*xn3  = 


roo 


' — 03 

r<x> 


Rf(x  + 0-  - e2)  h(e-| ) h(e2 )de -j d©2  (III-6) 

-ujECftt-epxJd  - r -e2)]h(01)h(02)de1de2  (III-7) 

-u2E[f(t  + r - 9-|)X^(t  - 92)]h(91)h(e2)d01d92  (III-8) 


where  Rf  is  the  autocorrelation  of  f(t)  and  h(e)  is  an  impulse-response 
function. 

To  make  the  perturbation  method  more  efficient,  we  have  to  simplify 
equations  (III-6),  (III-7),  and  (III-8).  Since  the  autocorrelation 
function  and  power  spectral  density  function  as  well  as  complex-frequency 
response  H(w)  and  unit-impulse  response  function  h(e)  are  related  through 
the  Fourier  Integral,  then  (III-6)  can  be  written  as 


E[X0X$]  = 


4>U)  | H (to)  |2 

« -00 


c()(ai)  | H (to)  |2  do) 


exp(io)t)  dw 


(IT  1-9) 
(II I -9a ) 


For  (III-7),  (II 1-8) , we  can  apply  the  Gaussian  property  [17]  and  get  the 
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final  result: 


E[X0XJ]  = E[X*X]  ] = -3co;0 


2 2 f <|)((d)  |H(ijo)  |2H‘ (ui)exp  (ion)  dm  (III-10) 

o 

When  u is  not  too  large,  for  simplicity,  we  evaluate  ( 1 1 1-3}  to  the 
first  order  of  . Substituting  ( 1 1 1 - 9 ) and  (III-10)  into  ( 1 1 1 -3) 

E[X  X*]  = 4> ( 03 ) | H(co)  |^exp  (iurr)  {1  - H'(u>)l  dw  (III-ll) 

n xQ 

(Ill-lla) 


2 _ 
ax  " 


where 


H(w)  = 


<j)(u>)  |H(a))  |2(1  - H 1 (o>) ) do 

o 0 


1 


2 2 + 9-Jr 

03^  - U)  + 


2 2 
ui  - w 


H'(u)  2 n2  2 

^ til  “ LiS  ) + 


y 


(2£u)nw)‘ 


Although  it  is  very  easy  using  the  digital  computer  to  calculate 
( I I I - 9a ) and  (Ill-lla),  the  pitfall  of  this  method  lies  in  the  fact 
that  we  cannot  prove  ( 1 1 1 -2)  is  convergent  and  actually  represents  a 
solution  of  (III-l)  within  an  acceptable  error.  Theoretically,  this 
method  is  only  for  slightly  nonlinear  systems,  i.e.  \x  is  very  small. 
However,  this  kind  of  knowledge  is  limited.  With  the  purpose  of 
broadening  the  applicability  of  this  method,  we  take  advantage  of  Monte 
Carlo  simulation  to  solve  (III-l)  and  compare  the  results  with  those 
from  the  perturbation  method. 


Monte  Carlo  Simulation 

We  have  already  known  that  the  excitation  f(t)  is  a stationary 
Gaussian  random  process  with  mean  zero  and  the  mean-square  spectral 
density  function  <j>(uj).  This  process  can  be  simulated  [18]  by  the  following 
series: 
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f0(t)  = | l cosUkt  + <j>k) 

k- 1 


(III-12) 


with  being  uniformly  distributed  between  0 and  2tt  and  independent 

of  4 . for  k f j and  with  <o.  being  a random  variable  independent  of 

a)-  for  k i j and  of  and  distributed  according  to  the  density  function: 

\i  J 


9(w) 


. U) 
2 

(T  f 


where 


°f  = 


4>((D)dn 


(III-13) 


(III-13a) 


is  the  variance  of  the  process  f(t). 

It  is  reasonable  to  simulate  f(t)  in  this  way  because  the  simulated 
process  fQ(t ) has  the  same  autocorrelation  function  and  spectral  density 
function  as  f(t). 

There  are  many  basic  methods  for  generating  variates  from  the  given 
probability  distribution.  The  inverse  transformation  method  seems  the 
most  convenient  if  it  can  be  employed.  We  first  obtain  the  cumulative 
distribution  function  FgU)  and  set  Fg(w)  - r: 


Fg(w)  = 


g(a>' ) doj'  = r 


(III-14) 


Now,  it  is  possible  to  find  a particular  u>.  corresponding  to  r • by  the 
inverse  function  of  F if  it  can  be  obtained 

= F'g  (ri)  (III-14a) 

where  r,  are  random  numbers  which  are  generated  by  different  processes 
[21].  In  case  we  cannot  express  in  terms  of  F~g(uO,  we  must  either 
calculate  a numberical  approximation  to  or  introduce  another  method 
instead  of  using  (III-14a). 

Once  the  generized  force  f(t)  is  simulated  by  (III-12)  and  (III-14a), 
(III-l)  can  be  solved  numerically.  The  accuracy  of  this  approach  is  checked 
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against  the  exact  solution  obtained  for  the  corresponding  linear  system. 

We  use  the  Duhamel  Integral  equation  [22]  to  solve  (III-l)  with  y = 0: 

X(t)  = A(t)sin(a)pt)  - B(t)cos(wpt)  (III-15) 

where 

*<*>  ■ if  Jo  V*'>  «ln-'5s) 

i ft  0Cwnt’ 

10  f0(,,)  JjfT  (m-i5b) 

Then,  we  obtain  the  variance  of  X(t)  and  compare  it  with  the  exact 
solution  from  (III-9a)  to  see  if  it  agrees  with  the  latter  in  an 
acceptable  range.  If  so,  we  may  conclude  the  simulation  is  good  and 
what  we  want  to  do  now  is  to  solve  (III-l)  with  y f 0,  i.e.  the  nonlinear 
system. 

Actually,  there  are  many  approximate  methods  to  solve  the  nonlinear 
equation  (III-l).  The  important  thing  is  that  we  have  to  avoid  large 
truncation  errors  which  result  in  a divergent  solution.  This  truncation 
error,  sometimes  called  the  discretization  error,  is  dependent  upon  the 
selection  of  the  interval  size.  We  must  choose  the  best  approach  to  meet 
both  economy  and  accuracy.  We  find  the  Runge-Kutta  formula  with  truncation 
error  of  fifth  order  to  be  a very  satisfactory  method  in  this  consideration. 
For  steps  smaller  than  some  critical  size  [23],  the  higher  error  order 
gives  less  error. 

Again,  we  compare  the  Runge-Kutta  solution  for  the  corresponding 
linear  system  with  that  from  (II 1-15).  When  these  studies  are  close, 
Runge-Kutta  may  be  applied  to  solve  (III-l).  The  variance  of  X(t)  will 
be  found. 
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Risk  Function  and  Reliability 

The  statistical  property  of  X(t)  has  been  obtained  by  the  above 
analysis.  Recalling  that  S(t)  = k(X(t)  + X3(t)),  we  obtain 


a2  = k2a2  + 6k2n/  + 15  k2p2a®  (III-16) 

S X X X 

With  this  information,  we  can  evaluate  the  reliability  of  the  structure 
if  we  know  the  initial  strength  characteristics  of  the  material. 

In  reliability  theory,  cessation  of  the  performance  of  mechanical 
functions  is  called  failure.  These  failures  can  be  separated  into  two 
groups.  First  is  the  failure  occurring  as  a result  of  stresses  exceeding 

their  limiting  or  critical  values.  The  failure  due  to  this  reason  is 

called  first-passage  failure.  Second  is  the  failure  due  to  repeated 
loadings,  i.e.  fatigue.  In  this  case,  the  loadings  are  not  large  enough 
to  cause  a failure  of  the  first  type  but  they  do  cause  the  successive 
incremental  reduction  of  the  limiting  response.  This  means  that 
repeated  loadings  hasten  the  reaching  of  first-passage.  If  the  particular 
structure  under  consideration  follows  the  power  law  for  flaw  propagation, 
(based  upon  the  Griffith-Irwin  equation  as  well  as  the  Extreme  Point  Process), 
[19]  we  get  the  risk  function,  the  probability  of  failure  after  a certain 
number  of  cycles,  is  as  follows: 


hp(n)  = 


“ exp(=4)  — ! 

0 2 as  /2tt*  an 


-(x-un)2 

exp  — 5 dx 

2 at 


(III-17) 


where 


2 _ 2 
an  ' aR. 


( 1 1 1 -1 7a ) 


2 

R„ 


being  the  mean  and  variance  of  the  initial 


(II 1-1 7b) 

resisting  stress  RQ 
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— 2 

respectively,  and  iu  and  cu  being  the  mean  and  variance  of  Zn, 
w n 3 n n 

Zn  = w E S.  with  S-  and  the  j-th  peak  of  S(t)  and  independent  of  S. 

n 2 j=i  j J k 

for  j f k also,  K is  a propagation  parameter. 

— 2 2 — 2 
Thus  far,  RQ,  are  known,  but  Uy  and  ay  are  unknowns.  To 

2 o n n 

find  Uy  and  a ..  , Crandall  [24]  first  considered  the  stress  process 

ln 

S(t)  to  consist  of  a sequence  of  half  cycles  each  of  duration  ir/u>n 
and  used  Rice's  definition  of  an  envelope: 

Uz  = J K [/2a.]3  r (1  + |)  (III~ 18a) 

z.n  c 

07  Tk  [^0  ]6  (1  - r2  (1  + |))  (II 1-1 8b) 

" M-l  2 

+ J-  K2[/?cs]6  r2  (1  + |)  x { 1 (M-k)[H(-f,-f;l,r*(T))-l]} 

K=  1 

Here,  H is  the  hypergeometric  function  and 

r* ( 1 ) = Ku  \ *M|H(M)|2cos(a,-u>n)Tda,)2  (III-19) 

as  0 

+ ( ^w)iH(.)|^in(_n)ldw)2]1/2 

^0  1 

4>(w)  |H(di)  is  the  spectral  density  of  the  process  X(t). 

p 

It  is  clear  that  the  further  evaluation  of  at  requires  the  value 

Ln 

of  r*(x)  to  specify  the  parameter  of  the  hypergeometric  function.  This 
is  a laborious  task.  One  way  to  overcome  this  difficulty  is  to  consider 
the  case  of  small  damping  and  to  assume  the  argument  of  the  hypergeometric 
function  is  oscillating  with  period  This  is  limited  to  white-noise 

type  excitations,  but  it  can  be  treated  as  an  approximate  value  for 
another  type  of  excitation,  especially  when  the  damping  is  small,  so 
the  complex- frequency  response  H(w)  is  very  sharply  peaked  at  about 
Thus  (II 1-1 8b)  is  approximated  as: 

a7  = 0.092  J K2[/2  a-]6  r2  (2.5) 
n * 5 


( 1 11-20) 
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Now  wo  are  able  to  calculate  (111-17).  However,  numerical 
integration  needs  thousands  of  operations  with  a high-speed  digital 
computer.  We  derive  the  following  formula  from  (1 11-17)  to  simplify 
the  evaluation: 


h (n)  = — ! ! 

p 2 Jl  anW 

1 , 1 


exp[-  \ (1  - -4-y)](/f  + erf  (—))  (III-21) 


2o 


2onW 


2°n 


where  W = 2a2  2 02  > and  erf  is  the  error  function  which  can  be  obtained 

by  the  asymptotic  series 


erf (x)  = 1 - -Le_x  (1  - -Ij  + ^ - ...) 
rV  x 2xJ  2 x 


(III-21a) 


Thus  the  reliability  of  a structure  which  will  survive  N stress  cycles 
is 


N-l 

UN)  = n [1  - h (n)3 
n=o  H 

N-l 

= exp[-  l h (n)]  for  h (n)  « 1 { II 1-22) 

n=o  H H 


A Numerical  Example  --  Wind  Pressure  on  a Structure 

As  an  example,  let  m = 1 kip  sec  /in  and  wn  = 1.732  rad/sec  for  the 
system  described  by  (1)  and  assume  the  system  starts  from  rest,  i.e.  with 
initial  conditions  X(t)  = X(t)  = 0.  The  system  is  excited  by  a strong 
gust  whose  spectrum  is  described  by  Davenport's  [25]  empirical  formula: 

Sy(f)df  = 4CV2 V473  dx  (II 1-23) 

( 1 "t  x ) 

where 

x = 4000  y (cycle/ foot) 

V = a standard  velocity  of  wind  at  33  feet  height 
C = a resistance  coefficient  of  the  ground  surface.  In  terrain 
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uniformly  covered  with  obstacles  30  - 50  feet  in  height, 
i.e.  residential  suburbs,  small  towns,  etc.  C - 0.015. 

So,  the  spectral  density  of  the  fluctuating  part  of  the  wind  pressure 
on  the  structure  may  be  assumed  as 


(pavcD)  2 

*(f)  = - - 4Cv 


(1  + x2)4'3 


(II 1-24) 


CD  = pressure  coefficient  =0.5  [26] 

Figure  6 shows  (II 1-24) . 

In  the  following  process  of  calculation,  we  must  always  recall  that 
w = 2irf , because  (II 1-24)  is  expressed  in  cycles,  but  in  the  previous 
analyses  the  formula  are  in  radians. 

From  (III-13) , (III-13a),  (III-14)  and  (III-14a),  we  get 


“i  = 2lrfi  = 2w  4000  [(1TF7)3  ' 1]1/2  1 “ l.-.-N  (II 1-25) 

With  the  aid  of  (III-12),  ( I I 1-1 3a ) , and  (II 1-25 ) , a typical  simulated 
excitation  with  wind  velocity  V = 120  ft/sec  and  N = 250  is  shown  in 
Fig.  7. 

1.  Monte  Carlo  Simulation 

Fig.  8 shows  the  comparison  between  different  approaches  for  the 
corresponding  linear  system.  The  triangles  indicate  the  exact  solution 
from  (III-9a)  and  crosses  denote  solutions  from  the  Duhamel  Equation 
(III-15)  based  upon  the  simulated  excitation  (III-12).  If  the  simulation 
(III-12)  is  good,  these  values  should  be  very  close.  Actually  the  technique 
of  simulation  (III-12)  is  dependent  on  the  selection  of  size,  the  number 
of  cosine  terms,  and  the  time  duration  for  computing  the  variance.  These 
factors  will  be  discussed  later.  From  Fig.  8,  we  may  conclude  that  the 
simulation  is  fairly  acceptable.  The  results  of  Fig.  8 at  the  same  time 
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indicate  the  accuracy  of  the  Runge-Kutta  Method  for  solving  ( I I 1-1 ) 
with  the  excitation  (III-12)  on  a step-by-step  timewise  basis.  Fig.  9 
shows  a typical  simulated  result  using  these  two  approximate  methods. 

The  fact  the  numerical  results  are  too  close  to  plot  as  two  distinct 
curves  reflects  the  feasibility  of  the  present  analysis.  The  computing 
time  for  this  approach  was  50.67  seconds. 

2.  Perturbation  Approach 

Figs.  10  to  15  show  a series  of  nonlinear  responses  each  with 
a different  value  of  y.  For  this  system:  (A)  If  y is  smaller  than  0.2 
the  perturbation  method  is  seen  to  be  reasonable  agreement  with  the 
Monte  Carlo  simulation  for  lower  values  of  excitations.  On  the  other 
hand,  some  simulated  solutions ■ (Figs.  10,  11,  and  12)  seem  to  be 
unreliable  because  their  values  are  higher  than  those  of  the  corresponding 
linear  system  shown  in  Fig.  3 and  reproduced  in  Figs.  10,  11,  and  12. 

This  is  due  to  the  fact  we  used  the  larger  time  interval  (0.4  sec.)  for 
higher  excitations  which  caused  a little  divergence  when  we  applied  the 
Runge-Kutta  formula.  This  deviation  can  be  reduced  when  the  nonlinearity 
increases.  (See  Figs.  13  and  14)  (B)  If  y is  larger  than  0.2,  application 
of  the  perturbation  method  will  cause  an  unacceptable  error,  as  shown  in 
Fig.  14.  However,  in  this  case,  the  simulated  approach  does  have  a great 
advantage.  All  data  are  plotted  in  Fig.  15. 

Fig.  16  describes  the  corresponding  linear  system  with  various 
damping  ratios.  It  is  quite  clear  that  the  higher  the  damping  ratio 
the  closer  these  two  methods  agree.  This  hints  that  the  perturbation 
method  will  do  well  when  y = 0.3  (or  larger)  and  simulated  solutions  will 
be  good  at  higher  excitations  even  with  a larger  interval  size  such  as 
0.4  sec.  These  are  showed  in  Figs.  17  and  18.  The  computing  time  for 
the  perturbation  approach  was  2.67  seconds. 
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Here,  considering  the  case  of  z - 0.02  and  v = 0.05  and  V = 100 
ft/sec,  (o2  = 0.5136  in2)  we  get  ox  = 0.284  in  from  Fig.  11.  Applying 
(III-16),  we  obtain: 

O2  = 0.744  in2  a.  = 0.862  in  (III-26) 

s s 

If  the  mean  of  the  initial  resistance,  RQ  is  3.45  kips  and  the  propagation 
parameter  K for  this  system  is  2.435  x 10  , with  the  aid  of  ( 1 1 1 - 1 8a ) , 

(II 1-20 ) , (III-21 ) and  (III-26),  we  obtain  the  estimates  of  the  risk 
function  h^( n ) which  is  the  curve  A in  Fig.  19.  Ine  reliability  is 
carried  out  by  ( 1 1 1-22 ) and  is  shown  by  curve  A in  Fig.  20.  Furthermore, 
we  consider  the  statistical  variation  of  the  initial  strength  RQ.  Curves 
B and  C in  Figs.  19  and  20  represent  the  results  for  variations  of  five 
and  ten  percent  respectively.  Similarly,  in  each  of  these  two  figures 
curve  D,  E and  F are  obtained  for  RQ  = 3.85  kips  and  variations  0,  five 
and  ten  percent  respectively. 

Conclusions 

Although  the  approaches  presented  are  developed  for  (III-l),  they  can 
be  extended  to  other  nonlinear  systems  by  the  same  proceudres.  From  the 
view  point  of  computing  for  the  specific  example  used,  the  perturbation 
method  is  more  efficient  than  the  Monte  Carlo  method.  The  perturbation 
method  takes  only  about  one-twentieth  the  computing  time  that  is  needed 
for  the  Monte  Carlo  simulation  technique.  In  general,  when  the 
perturbation  parameter  is  small,  the  method  works  well.  However,  if 
the  damping  ratio  increases  or  the  excitation  decreases,  the  method  still 
proves  to  be  adequate  even  if  the  perturbation  parameter  is  large.  This 
can  be  explained  by  the  fact  that  the  system  has  only  slight  nonlinearity. 
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Notjtjon_  used  in  Chapter  III 


X(t),  X,  X, 

i=l,Z...  1 

CL  ] 

s(t) 


sj 


n 

% 

T 

m 


Displacement  Response  Function 

Ensemble  average  or  expectation 
Stress  response  function 
The  j-th  peak  of  S(t) 

Natural  frequency  of  corresponding  linear  system 

Damped  vibration  frequency 

Fundamental  Period 

Mass  of  structure 

Damping  ratio 

Perturbation  of  nonlinear  parameter 


f(t) 

f0{t) 

<*(<*>) 


h(e) 

H(«) 

g(“) 

FgM 

F-]U) 


Stiffness  of  structure 
Initial  strength  of  structure 
Mean  of  Rq 

Standard  deviation  of  XQ,  X,  S(t),  and  R 
respectively  c 0 

Generalized  force,  excitation 
Simulation  of  f(t) 

Spectral  density  function  of  f(t) 

Variance  of  f(t) 

Autocorrelation  function  for  f(t) 
Impulse-response  function 
Complex- frequency  response 
Probability  density  function  for  w 
Cumulative  distribution  function  of  g(w) 
Inverse  function  of  FgU) 
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r,ri 


h (n) 
p 


H 

V 

L(N) 

Zn 

U^T 


C 

CD 

sv(f) 


Generating  random  numbers,  i ~ 1,2,...N 
Propagation  parameter 

Risk  function,  probability  of  structure's  failure  after 
n stress  cycles 

Hypergeometric  function 

Wind  velocity 

Reliability  of  structure 

n 3 
K/2  z S. 

Mean  of  ZR 

Variance  of  Zn 

Resistance  coefficient  of  ground  surface 
Pressure  coefficient 
Spectrum  of  wind  velocity 
Density  of  air,  0.0024  slug/ft3 
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CHAPTER  IV 

RESPONSE  OF  A NONLINEAR  SYSTEM  TO  NONSTATIONARY  RANDOM  EXCITATION 
Introduction 

The  transient  mean-square  response  of  a linear  single-degree-of- 
freedom  mechanical  system  to  certain  types  of  nonstationary  random 
excitation  has  been  studied  by  many  authors  [27,28,29,30].  The  nonstationary 
input  was  taken  in  the  form  of  a product  of  a well-defined  envelope 
function,  A(t)  and  a stationary  Gaussian  noise  with  zero  mean,  n(t). 

Caughey  and  Stumpf  [27]  have  examined  the  case  in  which  the  envelope 
function  A(t)  is  a unit  step  function  and  n(t)  is  assumed  to  be  either 
white  noise  or  broad-band  noise  whose  power  spectral  density  has  no 
sharp  peaks.  Results  of  their  analysis  were  applied  to  the  determination 
of  the  structural  response  to  earthquake  ground  motion.  Bolotin  [28] 
has  determined  the  mean-square  response  of  a structure  represented  by 
a second  order  differential  equation  to  earthquake  excitation.  In  his 
analysis,  he  considered  the  ground  acceleration  to  be  characterized  by 
the  product  of  an  exponentially  decaying  harmonic  correlation  function 
and  an  envelope  function,  A ( t ) = Ae 

In  a recent  paper  [29],  Barnoshi  and  Maurer  have  formulated  the  time 
varying  mean-square  response  of  a linear  si ngle-degree-of -freedom  system 
in  terms  of  the  system  frequancy  response  function  and  the  generalized 
spectral  density  function  of  the  input  excitation.  They  considered 
the  envelope  function  to  be  either  the  unit  step  function  or  a rectangular 
step  function.  Bucciarelli  and  Kuo  [30]  have  recently  obtained  an 
approximate  expression  for  the  mean-square  response  to  excitation 
characterized  by  a general  envelope  function  subject  only  to  the  restriction 
that  the  envelope  function  is  slowly  varying.  Their  work  also  gave  an 


est.im.iUxl  maximum  value  ot  the  moan-square  response. 

In  all  the  above  studies,  the  systems  treated  were  linear,  ihe 
present  study  presents  an  approximate  solution  to  the  problem  of 
transient  mean-square  response  of  a simple  nonlinear  system  to  a non- 
stationary random  excitation.  Only  systems  with  geometric  non! ineari ties 
(rather  than  materials  nonlinearities)  are  considered  and  the  nonlinear 
differential  equation  is  linearized  by  an  equivalent  linearization 
technique.  The  results  are  directly  applicable  to  determination  of 
response  of  an  elastic  flat  rectangular  or  circular  plate  subject  to 
random  lateral  loading  when  the  plate  is  approximated  as  a single- 
degree-of- freedom  system  characterized  by  its  central  deflection. 

Obviously  the  results  also  apply  to  other  one-degree-of-freedom  mechanical 
systems  excited  by  nonstationary  random  forces. 

Analysis 

Consider  a lightly  damped  single-degree-of- freedom  mechanical  system 
subjected  to  a random  excitation  and  governed  by  the  equation 

y(t)  + 2cainy(t)  + <^(y)  + g(y ) ) = f(t)  (IV-1) 

where 

C = fraction  of  critical  damping 

«n  = natural  frequency  of  the  corresponding  linear  system 
N 

g(y)  = z nky2k+i  uk  2 o (iv-2) 

The  nonstationary  random  excitation  f(t)  is  expressed  by 

f(t)  = A(t)n(t)  ( IV-3) 

where  A(t)  is  a well-defined  envelope  function  and  n(t)  is  a Gaussian 

stationary  random  process  with  zero  mean  and  autocorrelation  function  R (x). 

n' 
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2 

We  are  to  determine  the  mean-square  response  E[y  (t)]  to  an  input 
f(t)  when  the  envelope  function  is  a unit  step  function: 

A ( t ) = u(t)  (IV-4) 

and  n(t)  has  the  autocorrelation  functions 

Rn(  r)  = 2ttK06(t)  (IV-5) 

Although  various  methods  can  be  applied  to  determine  the  response 
of  nonlinear  systems,  the  equivalent  linearization  technique  will  be 
used  here.  This  technique  was  developed  by  Krylov  and  Bogoliuvov  for 
the  treatment  of  nonlinear  systems  under  deterministic  excitations, 
and  then  Booton  [31]  and  Caughey  [32]  applied  this  technique  to  problems 
of  random  vibrations. 

We  assume  that  an  approximate  solution  of  (IV-1)  can  be  obtained 
from  the  linearized  equation 

y + 2(iey  + *1  = f(t)  ( IV  -6) 

where  (3g  is  the  equivalent  linear  damping  coefficient  and  u>e  is  the 
equivalent  linear  stiffness.  The  error  "e"  due  to  linearization  is  given 
by  the  difference  between  (IV-1)  and  (IV-6),  i.e., 

e = 2(cwn  - ae)y  + (c^  - wj;)y  + g(y)c^  (IV-7) 

2 

The  variables  and  we  are  chosen  so  as  to  minimize  the  mean-square 

2 2*2' 
error  E[e  ].  The  resulting  values  involve  E[y  ],  [y  ] and  E[yy],  hence 

it  is  necessary  to  know  the  probability  density  function  p(y,y).  In 

general,  however,  p(y,y)  is  not  known.  If  the  input  is  Gaussian  and 

the  nonlinearities  of  the  system  are  small,  then  the  response  of  the 

linearized  equation  (II 1-6 ) is  also  Gaussian.  Therefore,  the  assumption 

is  made  that  the  probability  density  function  p(y,y)  is  Gaussian  with 

covariances  to  be  determined.  Before  constructing  the  probability  density 

function,  however,  it  is  necessary  to  find  ensemble  averages  of  y and  y 
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and  these  are  readily  found  by  DuhameVs  integral  to  be 

rt 


K[y]  - 


o 


h( t- I )t[ t( I ) Jd  t 


and  if  we  assume  that  E[f(t)]  = 0,  then 

E[y]  = o 

Similarly,  the  ensemble  average  of  y is  obtained  as 


E[y]  = j h(t-T)E[f(T)]dx 

= 0 


(1V-8) 


(IV-9) 


(IV-10) 


Thus,  the  assumed  Gaussian  probability  density  function  p(y,y)  takes 
the  form 

p(y,y)  ~TTo  exp(-ay2  + 2byy  - cy2)  (IV-11 ) 


2*(det(K))1/2 


where 

a = E[y2]/(2det(K)) 

b = E[yy]/(2det(K))  (IV-12) 

c = E[y2]/ (2det(K) ) 
det(K)  = E[y2]E[y2]  - (E[yy])2 
This  leads  to 


2l3e 


2*. 
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(E[y2])k 


( IV- 1 3 ) 
(IV-14) 


It  is  interesting  to  observe  that  the  above  equivalent  linear  damping 

2 

23e  and  stiffness  are  identical  to  those  found  for  a stationary  process 
in  which  case  E[yy]  is  equal  to  zero. 

If  the  nonlinearity  is  involved  only  in  the  velocity  term  such  as 
9 (y ) it. is  possible  to  demonstrate  thatthe  equivalent  linear  damping 
and  stiffness  for  a nonstationary  process  are  identical  to  those  for  a 
stationary  process. 
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In  all  cases,  the  mean-square  response  E[jr]  at  any  time  t is 

2 

obtained  from  the  expected  value  of  (y(t))  over  the  ensemble  response. 

Use  of  Duhamel's  integral  leads  to 

o i rt:  ft: 

E[y  ] " ~2  exp{-c^A)^(2t-T-T,)}sinuJd(t-T,)sinaJd(t-T)A(T)A('[,)  • 


ud  J0J0 


R ( T—  T1  JdldT1 


where 


N 


2 _ . 2/1  A (2k+1 )! 

2kk! 

If  the  input  is  white  noise,  then  (IV-15)  becomes 


4 = + i uk  (E[y2])k  - c2) 

d n k=1  2Kk! 


(IV-15) 


(IV-16) 
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— sin  (uji 

“d  d 

(IV-18) 


sin  wdt  + 77-^  sin  w^t}] 


Employing  (IV-16),  (IV-18)  becomes  a nonlinear  algebraic  equation  for 
E[y  ] since  ^ is  a function  of  E[y  ] in  (IV-16).  This  type  of  equation 
generally  has  more  than  one  solution.  However,  from  physical  considerations 
the  desired  solution  will  be  that  one  close  to  the  solution  of  the 
corresponding  linear  system  because  only  a weakly  nonlinear  system  is 
being  considered.  It  is  convenient  to  solve  this  system  through  use  of 
Newton* s method  of  tangents  together  with  iteration.  As  an  example,  let 
us  consider  the  case 

g(y)  = vy3  (iv-19) 

For  various  values  of  y and  damping  coefficient  c,  E[y2]  is  computed 
and  the  normalized  plots  such  as  shown  in  Figure  21  result.  The  normalization 
factor  is  determined  by  the  stationary  mean-square  response  of  the  linear 
system 


[54] 


E[yJ]  = "K0/4 


(IV-20) 


The  parameter  ji  is  chosen  in  such  a manner  that  given  n the  stationary 
mean- square  response  reaches  40  percent,  60  percent  and  80  percent 
Of  E[y2]$. 

If  the  damping  is  small,  (IV-18)  can  be  approximated  by 

E[y2]  * E[y2]  0,  - <L^1  (IV-21) 

0 S 1 + 3pE[yZ] 

from  which  the  following  approximate  solution  is  obtained: 

(IV-22) 

From  the  above  results  it  is  easy  to  demonstrate  that  the  transient 
mean-square  response  for  both  linear  and  nonlinear  systems  does  not  exceed 
the  stationary  mean-square  response  zo  white  noise. 

If  instead  of  (IV-5)  the  function  n(t)  has  the  autocorrelation 

( IV -23 ) 

(termed  correlated  noise)  then  (IV-15)  becomes 

K_e"2'>nt  rtft 

o‘o 


E[y2]  = jU  {[1  + 12ME[y2]s(l  - e_24‘»nt)]V2  . 1} 


Rn(i)  = K0  exp(-«|i|)cosp, 


E[y2]  - — 


exp[r>n(x+i  ' ) - |]A(r)A(t 1 )sina>d(t-i) 


sinwd(t-i' )cos&(x-t* )didi 1 (IV-24) 

Use  of  (IV-4)  in  (IV-24)  leads  to  the  desired  expression  for  E[y2].  In 
the  interest  of  brevity  this  lengthy  expression  is  not  given.  However, 
inspection  of  it  reveals  that  the  mean-square  response  depends  upon 
interrelationships  between  damping  c,  the  corresponding  linear  system 
natural  frequency  wn,  the  decay  constant  and  the  frequency  of  the 
correlation  function. 

For  a white  noise  input,  only  the  value  of  damping  of  the  system 
influences  the  time  required  to  attain  stationarity.  However,  for  the 
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correlated  noise  input,  the  time  required  for  the  response  to  reach 
a stationary  value  is  influenced  not  only  by  the  system  damping 
coefficient  c but  also  by  the  decay  constant  « of  the  input  noise. 

Inspection  of  results  indicates  that  as  a decreases,  i.e.  the  power 
spectral  density  has  a sharp  peak  at  some  frequency,  then  the  transient 
response  tends  to  exceed  the  stationary  value.  Another  interesting 
result  is  that  the  nonlinear  response  becomes  greater  than  the  corresponding 
linear  response  under  certain  conditions  even  if  the  system  has  hardening 
spring-type  nonlinearity.  An  example  of  this  is  shown  in  Figure  22. 

Conclusions 

The  time  varying  mean-square  response  of  a nonlinear  single-degree- 
of-freedom  mechanical  system  to  nonstationary  random  excitation  characterized 
by  the  product  of  an  envelope  function  and  a stationary  Gaussian  random 
process  has  been  considered.  A unit  step  envelope  function  was  considered 
in  conjunction  with  both  correlated  and  white  noise  with  zero  mean.  The 
nonlinear  governing  equation  was  linearized  by  the  method  of  equivalent 
linearization. 

For  the  nonstationary  process  it  has  been  shown  that  the  equivalent 
linear  damping  coefficient  and  the  equivalent  linear  stiffness  for  the 
system  with  nonlinearities  involved  only  in  displacements  or  only  in 
velocities  are  the  same  as  those  for  the  stationary  process. 

The  mean-square  response  depends  upon  the  coefficients  of  the  system 
equation,  the  shape  of  the  envelope  function,  and  the  parameters  of  the 
autocorrelation  of  the  process  n(t).  It  was  proved  that  for  white  noise 
modulated  by  a unit  step  function,  the  transient  mean-square  response  never 
exceeds  the  stationary  response.  However,  the  mean-square  response  to 


corrected  noise  modulated  by  a unit  stop  function  may  exceed  its 
stationary  value,  especially  when  the  power  spectral  density  of  the 
process  n(t)  has  a sharp  peak,  and  its  maximum  value  becomes  several 
times  the  stationary  value. 

It  has  also  been  shown  that  the  mean-square  response  of  the  system 
with  cubic  hardening  spring-type  nonlinearity  may  be  greater  than  the 
corresponding  linear  system  response  under  certain  conditions. 

The  results  presented  are  directly  applicable  to  determination  of 
response  of  a single-degree-of-freedom  mechanical  system  to  nonstationary 
random  excitation.  In  the  case  of  a flat  rectangular  plate  subject  to 
random  lateral  loading  the  response  at  the  center  of  the  plate  is 
frequently  the  one  of  greatest  interest  and  in  this  case  the  vibrating 
plate  may  be  represented  as  a one-degree-of-freedom  system. 
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Fig.  6.  Spectral  density  function  of  wind  pressure 


[61] 


RMS  response  of  system  to  excitation  withfs  0.02  and 
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Pig.  9.  A selection  of  typical  simulated  response 


Simulation  (i£q.  12 


Fig.  10.  RMS  response  of  system  to  excitation  with  £=  0.02 


Simulation  (£q.  12) 


response  of  system  to  excitation  with  £s  0.02  and 


Fig.  12,  RMS  response  of  system  to  excitation  with  £ = 0.02  an 


response  of  system  to  excitation  with  t = 0.02 


Exact  solution  (Eq.  11a) 


Fig.  16.  rms  response  of  system  to  excitation  with  u = 0 (linear  system) 


Perturbation  (2q.  lib) 
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response  of  system  to  excitation  with  £ = 0.05 


Perturbation  (2q.  11b) 
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Fig.  18.  RMS  response  of  system  to  excitation  with  % = 0.10  and  % =.  0.20 
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Risk  function  as  a function  of  cycles 


1,0 
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Reliability  of  structure  as  a function  of  cycles 


Normalized  Mean-square  Response  E(y  ]/K  C 
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Mean-square  Response  of  th*  Nonlinear  Systems  with  verlpus  Nonl inear it i*e 
to  White  Noise  Modulated  by  a Unit  Step  Function.  Systpm  Damping  C f 0.0S 

Figure  21 


Mean-square  Response  of  the  Systems  to  Correlated  Noise  modulated  by  a Unit  Step  function. 
System  Damping  C ■ 0.1 


Figure  22 


